library(yardstick)
library(DALEXtra)
library(pdp)
library(themis)
library(ranger)
library(dplyr)
library(ggplot2)
library(kknn)
library(naivebayes)
library(glmnet)
library(discrim)
library(tidyverse)
library(purrr)
library(tidymodels)
library(visdat)
library(kernlab)
library(corrplot)
library(baguette)

# data <- read_csv("aprendizagem_maquina/diabetes.csv") |> mutate(Outcome = as.factor(Outcome))
data <- read_csv("diabetes.csv") |> mutate(Outcome = as.factor(Outcome))
cols <- c("Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI")
data[cols][data[cols] == 0] <- NA


Ausentes <- data |>
  vis_dat() +
  labs(
    y = "Observações"
    ) +
  scale_fill_manual(
    values = viridis::mako(10)[5:6],
    labels = c("Qualitativa", "Quantitativa", "Ausentes")
    ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
  )

set.seed(123)
data_split <- initial_split(data, prop = 0.8, strata = Outcome)
train_data <- training(data_split)
test_data <- testing(data_split)

rec <- recipe(Outcome ~ ., data = train_data) %>%
  step_impute_median(all_numeric_predictors())

prep_rec <- prep(rec)
train_data_eda <- bake(prep_rec, new_data = train_data)
test_data_eda <- bake(prep_rec, new_data = test_data)

Frequência <- train_data_eda |>
  ggplot(aes(x = forcats::fct_infreq(as.factor(Outcome)))) +
  geom_bar(
        aes(
          y = after_stat(count) / sum(after_stat(count)),
          fill = Outcome
          ),
        show.legend = FALSE,
  ) +
  theme_bw() +
  theme(axis.title.x = element_blank()) +
  scale_fill_manual(
    values = viridis::mako(5)[3:5]
    ) +
  labs(y = NULL, x = NULL)

Distribuição <- train_data_eda |>
  pivot_longer(
    !Outcome,
    names_to = "variable",
    values_to = "values"
  ) |>
  mutate(Outcome = as.factor(Outcome)) |>
  ggplot(
    aes(
      x = values,
      fill = Outcome
      )
    ) +
  geom_violin(
    aes(y = variable),
    draw_quantiles = c(0.25, 0.5, 0.75),
    ) +
  labs(x = "", y = "") +
  facet_wrap(vars(variable), scales = "free") +
  scale_fill_manual(
    values = viridis::mako(5)[3:5],
    name = "Outcome"
    ) +
  theme_bw() +
  theme(axis.text.y = element_blank())

data_cor <- train_data_eda |>
  (\(x) {
     scaled <- x |>
       select_if(is.numeric) |>
       scale() |>
       as_tibble()

     x |>
       select_if(is.character) |>
       cbind(scaled)
  })()

data_cor |>
  select_if(is.numeric) |>
  cor() |>
  corrplot(
      method="circle",
      col = viridis::mako(200),
      diag = FALSE,
      addgrid.col = "#00000040",
      tl.pos = "l",
      tl.cex = 0.6,
      tl.col = "black",
      tl.offset = 0.25,
      cl.length = 9,
      cl.cex = 0.6,
      cl.ratio = 0.25,
      family = "mono"
  )

set.seed(123)

cv_folds <- vfold_cv(train_data,
                     v = 10,
                     strata = Outcome
                     )

knn_spec <- nearest_neighbor(neighbors = tune(), weight_func = tune()) |>
  set_engine("kknn") |>
  set_mode("classification")

nbayes_spec <- naive_Bayes(smoothness = tune(), Laplace = tune()) |>
  set_engine("naivebayes") |>
  set_mode("classification")

lda_spec <- discrim_linear() |>
  set_engine("MASS") |>
  set_mode("classification")

qda_spec <- discrim_quad() |>
  set_engine("MASS") |>
  set_mode("classification")

logistic_spec <- logistic_reg(
  engine = "glmnet",
  penalty = tune(),
  mixture = tune()
)

dt_spec <- decision_tree(
  cost_complexity = tune(),
  min_n = tune(),
  tree_depth = tune()) |>
  set_engine(engine = "rpart") |>
  set_mode("classification")

rf_spec <- rand_forest(mtry = tune(),
                       min_n = tune(),
                       trees = tune()) |>
  set_engine(engine = "ranger") |>
  set_mode("classification")

bt_spec <- bag_tree(
        cost_complexity = tune(),
        min_n = tune(),
        tree_depth = tune()
    ) |>
    set_engine("rpart") |>
    set_mode("classification")

rf_spec <- rand_forest(
        mtry = tune(),
        min_n = tune(),
        trees = tune()
    ) |>
    set_engine(engine = "ranger") |>
    set_mode("classification")

xgb_spec <- boost_tree(
    tree_depth = tune(),
    learn_rate = tune(),
    loss_reduction = tune(),
    min_n = tune(),
    sample_size = tune(),
    trees = tune(),
    mtry = tune()
    ) |>
    set_engine(engine = "xgboost") |>
    set_mode("classification")

lsvm_spec <- svm_linear(cost = tune(),
                        margin = tune()) |>
  set_engine(engine = "kernlab") |>
  set_mode("classification")

rsvm_spec <- svm_rbf(cost = tune(),
                     rbf_sigma = tune(),
                     margin = tune()) |>
  set_engine(engine = "kernlab") |>
  set_mode("classification")

psvm_spec <- svm_poly(cost = tune(),
                      degree = tune(),
                      scale_factor = tune(),
                      margin = tune()) |>
  set_engine("kernlab")  |>
  set_mode("classification")

mlp_spec <- mlp(
  hidden_units = tune(),
  penalty = tune(),
  epochs = tune()
) |>
  set_engine("nnet") |>
  set_mode("classification")

wf = workflow_set(
  preproc = list(rec),
  models = list(
        knn_fit = knn_spec,
    nbayes_fit = nbayes_spec,
    lda_fit = lda_spec,
    qda_fit = qda_spec,
    linear_fit = logistic_spec,
        dt_fit = dt_spec,
        bt_fit = bt_spec,
        rf_fit = rf_spec,
        xgb_fit = xgb_spec,
    lsvm_fit = lsvm_spec,
    rsvm_fit = rsvm_spec,
    psvm_fit = psvm_spec,
    mlp_fit = mlp_spec
    )
  )  |>
  mutate(wflow_id = gsub("(recipe_)", "", wflow_id))


grid_ctrl = control_grid(
  save_pred = TRUE,
  parallel_over = "resamples",
  save_workflow = TRUE
)

grid_results = wf  |>
  workflow_map(
    seed = 123,
    resamples = cv_folds,
    grid = 50,
    control = grid_ctrl
  )

autoplot(grid_results, metric = "roc_auc") +
  theme_bw() +
  labs(y = "Métrica", x = "")

autoplot(grid_results, select_best = TRUE) +
  theme_bw() +
  labs(y = "Métrica", x = "")

results_acc = workflowsets::rank_results(
  grid_results,
  select_best = TRUE,
  rank_metric = "roc_auc"
  ) |>
  filter(.metric == "roc_auc") |>
  dplyr::select(wflow_id, mean, std_err, model, rank)

model_ids <- c(
  "linear_fit", "knn_fit",
  "nbayes_fit", "lda_fit",
  "qda_fit", "dt_fit",
  "bt_fit", "rf_fit", "xgb_fit",
  "rsvm_fit", "lsvm_fit", "psvm_fit", "mlp_fit")

best_sets <- map(model_ids, ~ grid_results %>%
                   extract_workflow_set_result(.x) %>%
                   select_best(metric = "roc_auc"))

names(best_sets) <- model_ids

test_results <- function(rc_rslts, fit_obj, par_set, split_obj) {
  rc_rslts %>%
    extract_workflow(fit_obj) %>%
    finalize_workflow(par_set) %>%
    last_fit(
      split = split_obj,
      metrics = metric_set(accuracy, roc_auc, brier_class
      )
    )
}

test_results_list <- map2(model_ids, best_sets,
                          ~ test_results(grid_results, .x, .y, data_split))

metrics_table <- map_dfr(test_results_list, collect_metrics, .id = "model_id") |>
  dplyr::select(model_id, .metric, .estimate) |>
  pivot_wider(names_from = .metric, values_from = .estimate) |>
  mutate(across(where(is.numeric), round, 4)) |>
  mutate(
    Modelo = c("Logistic", "K-NN",
    "Naive Bayes", "LDA",
    "QDA", "Árvore de Decisão", "Boost tree",
    "Floresta aleatória", "XGBoost",
    "RSVM", "LSVM", "PSVM", "Neural Network"
    )
  ) |>
  relocate(Modelo) |>
  dplyr::select(-model_id)

colnames(metrics_table) <- c("Modelo", "Accuracy", "ROC AUC", "Brier Score")

metrics_table <- as_tibble(metrics_table)

metrics_table |>
  arrange(desc(`ROC AUC`)) |>
  knitr::kable()

rec <- recipe(Outcome ~ ., data = train_data) |>
  step_impute_median(all_numeric_predictors()) |>
  step_mutate(across(all_numeric_predictors(), log1p)) |>
  step_mutate(
    N0 = BMI * SkinThickness,
    N8 = Age / Insulin,
    N13 = Glucose / DiabetesPedigreeFunction,
  ) |>
  step_smote(Outcome) |>
  step_nzv(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())

set.seed(123)

cv_folds <- vfold_cv(train_data,
                     v = 10,
                     strata = Outcome
                     )

knn_spec <- nearest_neighbor(neighbors = tune(), weight_func = tune()) |>
  set_engine("kknn") |>
  set_mode("classification")

nbayes_spec <- naive_Bayes(smoothness = tune(), Laplace = tune()) |>
  set_engine("naivebayes") |>
  set_mode("classification")

lda_spec <- discrim_linear() |>
  set_engine("MASS") |>
  set_mode("classification")

qda_spec <- discrim_quad() |>
  set_engine("MASS") |>
  set_mode("classification")

logistic_spec <- logistic_reg(
  engine = "glmnet",
  penalty = tune(),
  mixture = tune()
)

dt_spec <- decision_tree(
  cost_complexity = tune(),
  min_n = tune(),
  tree_depth = tune()) |>
  set_engine(engine = "rpart") |>
  set_mode("classification")

rf_spec <- rand_forest(mtry = tune(),
                       min_n = tune(),
                       trees = tune()) |>
  set_engine(engine = "ranger") |>
  set_mode("classification")

bt_spec <- bag_tree(
        cost_complexity = tune(),
        min_n = tune(),
        tree_depth = tune()
    ) |>
    set_engine("rpart") |>
    set_mode("classification")

rf_spec <- rand_forest(
        mtry = tune(),
        min_n = tune(),
        trees = tune()
    ) |>
    set_engine(engine = "ranger") |>
    set_mode("classification")

xgb_spec <- boost_tree(
    tree_depth = tune(),
    learn_rate = tune(),
    loss_reduction = tune(),
    min_n = tune(),
    sample_size = tune(),
    trees = tune(),
    mtry = tune()
    ) |>
    set_engine(engine = "xgboost") |>
    set_mode("classification")

lsvm_spec <- svm_linear(cost = tune(),
                        margin = tune()) |>
  set_engine(engine = "kernlab") |>
  set_mode("classification")

rsvm_spec <- svm_rbf(cost = tune(),
                     rbf_sigma = tune(),
                     margin = tune()) |>
  set_engine(engine = "kernlab") |>
  set_mode("classification")

psvm_spec <- svm_poly(cost = tune(),
                      degree = tune(),
                      scale_factor = tune(),
                      margin = tune()) |>
  set_engine("kernlab")  |>
  set_mode("classification")

mlp_spec <- mlp(
  hidden_units = tune(),
  penalty = tune(),
  epochs = tune()
) |>
  set_engine("nnet") |>
  set_mode("classification")

wf = workflow_set(
  preproc = list(rec),
  models = list(
        knn_fit = knn_spec,
    nbayes_fit = nbayes_spec,
    lda_fit = lda_spec,
    qda_fit = qda_spec,
    linear_fit = logistic_spec,
        dt_fit = dt_spec,
        bt_fit = bt_spec,
        rf_fit = rf_spec,
        xgb_fit = xgb_spec,
    lsvm_fit = lsvm_spec,
    rsvm_fit = rsvm_spec,
    psvm_fit = psvm_spec,
    mlp_fit = mlp_spec
    )
  )  |>
  mutate(wflow_id = gsub("(recipe_)", "", wflow_id))


grid_ctrl = control_grid(
  save_pred = TRUE,
  parallel_over = "resamples",
  save_workflow = TRUE
)

grid_results = wf  |>
  workflow_map(
    seed = 123,
    resamples = cv_folds,
    grid = 50,
    control = grid_ctrl
  )

results_acc = workflowsets::rank_results(
  grid_results,
  select_best = TRUE,
  rank_metric = "roc_auc"
  ) |>
  filter(.metric == "roc_auc") |>
  dplyr::select(wflow_id, mean, std_err, model, rank)

model_ids <- c(
  "linear_fit", "knn_fit",
  "nbayes_fit", "lda_fit",
  "qda_fit", "dt_fit",
  "bt_fit", "rf_fit", "xgb_fit",
  "rsvm_fit", "lsvm_fit", "psvm_fit", "mlp_fit")

best_sets <- map(model_ids, ~ grid_results %>%
                   extract_workflow_set_result(.x) %>%
                   select_best(metric = "roc_auc"))

names(best_sets) <- model_ids

test_results <- function(rc_rslts, fit_obj, par_set, split_obj) {
  rc_rslts %>%
    extract_workflow(fit_obj) %>%
    finalize_workflow(par_set) %>%
    last_fit(
      split = split_obj,
      metrics = metric_set(accuracy, roc_auc, brier_class
      )
    )
}

test_results_list <- map2(model_ids, best_sets,
                          ~ test_results(grid_results, .x, .y, data_split))

metrics_table <- map_dfr(test_results_list, collect_metrics, .id = "model_id") |>
  dplyr::select(model_id, .metric, .estimate) |>
  pivot_wider(names_from = .metric, values_from = .estimate) |>
  mutate(across(where(is.numeric), round, 4)) |>
  mutate(
    Modelo = c("Logistic", "K-NN",
    "Naive Bayes", "LDA",
    "QDA", "Árvore de Decisão", "Boost tree",
    "Floresta aleatória", "XGBoost",
    "RSVM", "LSVM", "PSVM", "Neural Network"
    )
  ) |>
  relocate(Modelo) |>
  dplyr::select(-model_id)

colnames(metrics_table) <- c("Modelo", "Accuracy", "ROC AUC", "Brier Score")

metrics_table <- as_tibble(metrics_table)

metrics_table |>
  arrange(desc(`ROC AUC`)) |>
  knitr::kable()